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Variational Multiscale Stabilization and the 
Exponential Decay of Fine-scale Correctors 


Daniel Peterseim* 


Abstract This paper reviews the variational multiscale stabilization of standard fi¬ 
nite element methods for linear partial differential equations that exhibit multiscale 
features. The stabilization is of Petrov-Galerkin type with a standard finite element 
trial space and a problem-dependent test space based on pre-computed fine-scale 
correctors. The exponential decay of these correctors and their localisation to local 
cell problems is rigorously justified. The stabilization eliminates scale-dependent 
pre-asymptotic effects as they appear for standard finite element discretizations of 
highly oscillatory problems, e.g., the poor approximation in homogenization 
problems or the pollution effect in high-frequency acoustic scattering. 


1 Introduction 

In the past decades, the numerical analysis of partial differential equations (PDEs) 
was merely focused on the numerical approximation of sufficiently smooth solutions 
in the asymptotic regime of convergence. In the context of multiscale problems (and 
beyond), such results have only limited impact because the numerical approximation 
will hardly ever reach the asymptotic idealised regime under realistic conditions. Al¬ 
though a method performs well for sufficiently fine meshes it may fail completely on 
coarser (and feasible) scales of discretization. This happens for instance if the PDE 
exhibits rough and highly oscillatory solutions. Among the prominent applications 
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are the numerical homogenization of elliptic boundary value problems with highly 
varying non-smooth diffusion coefficient, high-frequency time-harmonic acoustic 
wave propagation, and singularly perturbed problems such as convection-dominated 
flow. 

The numerical approximation of such problems by finite element methods (FEMs) 
or related schemes is by no means straight-forward. The pure approximation (e.g. 
interpolation) of the unknown solutions by finite elements already requires high spa¬ 
tial resolution to capture fast oscillations and heterogeneities on microscopic scales. 
When the function is described only implicitly as the solution of some partial dif¬ 
ferential equation, its approximation faces further scale-dependent pre-asymptotic 
effects caused by the under-resolution of relevant microscopic data. Examples are 
the poor approximation in homogenization problems (see Fig. 1) and the pollu¬ 
tion effect [BSOO] for Helmholtz problems with large wave numbers (see Fig. 2). We 
shall emphasise that, in the latter case, the existence and uniqueness of numerical 
approximations may not even be guaranteed in pre-asymptotic regimes. 

Such situations require the stabilization of standard methods so that eventually 
a meaningful approximation on reasonably coarse scales of discretization becomes 
feasible. This paper aims to present a general framework for the stabilization of 
FEMs for multiscale problems with the aim to significantly reduce or even elimi¬ 
nate pre-asymptotic effects due to under-resolution. Our starting point will be the 
Variational Multiscale Method (VMS) originally introduced in [Hug95, HFMQ98]. 
The method provides an abstract framework how to incorporate missing fine-scale 
effects into numerical problems governing coarse-scale behaviour [HS07]. One may 
interpret the VMS as a Petrov-Galerkin method using standard EE trial spaces and 
an operator-dependent test space that needs to be precomputed in general. 



Fig. 1 Failure of FEM in homogenization problems: Consider the periodic problem 
(x) ^ U£ (x) = 1 in the unit interval with homogeneous Dirichlet boundary condition, where 
Ae{x) := (2 + cos(27rx/e))“^ for some small parameter e > 0. The solution Ue = 4{x — x^) — 
4e sin(27r|) — ^xsin(27r|) — ^ cos(27r|) + is depicted in blue for e = 2~^. The FI¬ 
FE approximation (o) on a uniform mesh of width h interpolates the curve x ^ 2a/3(x — x^) when¬ 
ever h is some multiple of the characteristic length scale e and, hence, fails to approximate Ue in 
any reasonable norm in the regime h>£. 
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The construction of this operator-dependent test space is based on some stable 
projection onto the standard FE trial space and a corresponding scale decomposi¬ 
tion of a function into its FE part given by the projection (the macroscopic/coarse- 
scale part) and a remainder that lies in the kernel of the projection operator (the 
microscopic/fine-scale part). The test functions are computed via a problem-depen- 
dent projection of the trial space into the space of fine-scale functions. This requires 
the solution of variational problems in the kernel of the projection - the fine-scale 
corrector problems. It has been observed empirically in certain applications that the 
Green’s function associated with these fine-scale corrector problems - the so-called 
fine-scale Green’s function [HFMQ98] - may exhibit favourable exponential de¬ 
cay properties [Mal05, HFMQ98] even though the decay of the classical full scale 
Green’s function is only algebraic. It is this exponential decay property that allows 
one to turn the VMS into a feasible numerical method [Mal05, LM07]. 

Very recently, the exponential decay was rigorously proved for the first time in 
[MP14] in the context of multi-dimensional numerical homogenization. A key in¬ 
gredient of the proof of [MP14] is the use of a (local) quasi-interpolation operator 
for the scale decomposition. Although the method of [MP14] still fits into the gen¬ 
eral framework of the VMS, it uses a different point of view on the method based 
on the orthogonalisation of coarse and fine scales with respect to the inner product 
associated with a symmetric and coercive model problem. This is why the method is 
now referred to as the Localized Orthogonal Decomposition (LOD) method. Subse¬ 
quent work showed that the ideas of [MP14] can be generalised to other discretiza¬ 
tion techniques such as discontinuous Galerkin [EGM13, EGMP13, Elf 15], Petrov- 
Galerkin formulations [EGH15], mixed methods [HHM15] and mesh-free methods 
[HMP15]. Moreover, the method can also be reinterpreted in terms of the multi¬ 
scale finite element method with special oversampling [HP 13]. The class of prob¬ 
lems that have been analysed by now includes semi-linear problems [HMP14b], 
high-contrast problems [PS 15, BP14], rough boundary conditions [HM14], prob¬ 
lems on complicated geometries [ELM 15], linear and non-linear eigenvalue prob- 



Ux = exp(—/o) is depicted in blue for K = 2^. The Pl-FE approximation (o) on a uniform mesh of 
width h = 2“^ > 6 • (wave length) fails to approximate due to the accumulation of phase errors. 
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lems [MP15b, HMP14a, MP15c], parabolic problems [MP15a], wave propagation 
[AH14, Petl4, GP15] and parametric problems [AH15]. 

This survey aims to reinterpret all those results in the abstract stabilization frame¬ 
work of the original VMS (Section 2) and aims to illustrate how the exponential 
decay of the fine-scale Green’s function can be quantified (Section 3). We will 
show how these abstract results lead to super-localised numerical homogenization 
[MP14, HP13] (Section 4) and pollution-free time-harmonic acoustic scattering 
(Section 5) [Petl4, GP15]. Section 6 contains some final remarks and also also 
identifies methodological similarities and differences with some other numerical ap¬ 
proaches that receive great attention these days, e.g., discontinuous Petrov-Galerkin 
methods (dPG) [DGll], Trefftz-type methods [GHP09] and Isogeometric Analysis 
(IGA) [CHB09]. 


2 Abstract variational multiscale stabilization 


This section is concerned with an abstract variational problem in a complex Hilbert 
space V as it appears for the weak formulation of second order PDEs. In this context, 

V is typically some closed subspace of the Sobolev space for some 

bounded Lipschitz domain Q Let a denote a bounded sesquilinear form on 

V xV and let F eV' denote a bounded linear functional on V. We wish to find ueV 
satisfying the linear variational problem 

\/v eV \ a{u^v) = F{v). (1) 


We assume that the sesquilinear form a satisfies the inf-sup condition 


a := inf sup 


a{v,w) 


= inf sup 


a{v,w) 


lrl|v||w||v lrl|v||w||v 


> 0 . 


( 2 ) 


Under this condition, the abstract problem (1) is well-posed, i.e., for all L G V' there 
exists a unique solution u eV and the a priori bound 


holds true; see, e.g., [Bab71]. 

We wish to approximate the unknown solution u of (1) by some computable 
function. The standard procedure for approximation is the Galerkin method which 
simply chooses a finite-dimensional subspace Vh CV (that contains simple func¬ 
tions such as piecewise polynomials) and restricts the variational problem (1) to 
this subspace. Usually, Vh belongs to some family of spaces parametrised by some 
abstract discretization parameter H, for instance the mesh size. This parameter (or 
set of parameters) provides some control on the approximation properties of Vh as 
^ ^ 0 at the price of an increasing computational cost in the sense of dimVn 
The Galerkin method seeks a function Ghu G Vh satisfying 
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^vh ^Vh ' ci{Ghu,vh) = F{vh) { = a{u,VH))- (3) 


Recall that the well-posedness of the original problem (1) does not imply the well- 
posedness of the discrete variational problem (3) but needs to be checked for the 
particular application via discrete versions of the inf-sup condition (2). In many 
cases, such conditions are only satisfied for H sufficiently small. This means that 
there is some threshold complexity for computing any Galerkin approximation and 
this threshold can be out of reach. Even if a Galerkin solution Gru exists and is 
computable, it might not provide the desired accuracy or does not refiect the relevant 
characteristic features of the solution, as we have seen in the introduction. 

Therefore, we are interested in computing projections onto the discrete space Vr 
other than the Galerkin projection. Let Ir :V ^Vr denote such a linear surjective 
projection operator and let us assume that it is bounded in the sense of the space 
^(V) of linear operators from V to V with finite operator norm 


¥h\\^{V) ■= sup 


llfevllv 

l|v|k 


oo. 


Implicitly, we also assume that this operator norm does not depend on the discretiza¬ 
tion parameter // in a critical way. Possible choices of 1r include the orthogonal 
projection onto Vr with respect to the inner product of V or any Hilbert spaces 
L D y containing V and mainly (local) quasi-interpolation operators of Clement or 
Scott-Zhang type as they are well-established in the finite element community in the 
context of a posteriori error estimation [Cle75, SZ90, Car99, CV99]. 


2.1 Petrov-Galerkin characterisation of finite element projections 

The Galerkin projection Gr is designed in such a way that its computation requires 
only the known data F associated with the unknown solution u. This section charac¬ 
terises the projection /// G ^(V) as a Petrov-Galerkin discretization of (1) using Vr 
as the trial space and a non-standard test space Wr C V that depends on the problem 
and the projection. The definition of Wr rests on the trivial observation that, for any 

vGV, _ 

a{lRU,v) = F{v) — a{u — Iru^v). (4) 

The choice of a test function v in the subspace 

Wr :={weV I Vz G KctIr : a{z,w) = 0} (5) 

annihilates the second term on the right-hand side of (4) and, hence, 

ci{Iru,wr) = F{wr) 

holds for all wr G Wr. This shows that Iru is a. solution of the Petrov-Galerkin 
method: Find ur G Vr such that 



6 


Daniel Peterseim 


\/wh G Wh : ci{uh,wh) = F{wh)- ( 6 ) 

This characterisation of Ih is well known from the variational multiscale method as 
it is presented in [HS07]. 

The question whether or not (6) has a unique solution can not be answered under 
the general assumptions made so far. We need to assume the missing uniqueness to 
be able to proceed and one way of doing this is to assume that the dimensions of 
trial and test space are equal, 


dim Wh = dim Vh- (7) 

In the present setting with a bounded operator Ih, this condition is equivalent to 
the well-posedness of the discrete variational problem (6), i.e., it admits a unique 
solution uh = Ihu ^Vh and 


\\uh\\v < ||^//||^( l /) l | M|k < 

The a priori estimate in turn implies a lower bound of the discrete inf-sup constant 
of the Petrov-Galerkin method by the quotient of the continuous inf-sup constant a 
and the continuity constant of Ih, 

. a{vH,WH) . a ^ ■ f a{vH,WH) 

\\^h\\v\\^^h\\v \Vh\\^(v) 11^//IIV Ik//IIV 


The test space Wh is the ideal test space for our purposes in the following sense. 
Assuming that we have access to it, the method (6) would enable us to compute Ihu 
without the explicit knowledge of u. Although this will rarely be the case, we will 
see later that Wh can be approximated very efficiently in relevant cases. The discrete 
inf-sup conditions then indicate that the sufficiently accurate approximation of Wh 
will not harm the method, its stability properties or its subsequent error minimisation 
properties. 

The continuity of the projection operator Ih readily implies the quasi-optimality 
of the Petrov-Galerkin method (6), 

11^ —^//||y = 11(1 —7//)w||v < ||///|k(y) min \\u — vh\\v‘ (8) 

vh^Vh 

Here, we have used that ||///|k(v) = ||1 ~f//|k(v); see e.g. [Szy06]. More impor¬ 
tantly, the same arguments show that the Petrov-Galerkin method is quasi-optimal 
with respect to any other Hilbert space L D V with norm \\'\\l whenever Ih G ^(L), 

\\u-uh\\l<\\Ih\\^{l) min \\u-vh\\l‘ 

vh^Vh 

These quasi-optimality make the ansatz very appealing and motivates further in¬ 
vestigation. Hence, in the remaining part of the paper, it is our aim to turn the method 
into a feasible numerical scheme while preserving these properties to a large extent. 
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Although the discrete stability of the method depends on the stability properties of 
the original problem and, hence, on parameters such as the frequency in scatter¬ 
ing problems, the quasi-optimality depends only on 1h and not necessarily on the 
problem. 


2.2 Characterisation of the ideal test space 


A practical realisation of the Petrov-Galerkin method (6) requires a choice of bases 
in the discrete trial Vh and test space Wh- These choices will have big impact on the 
computational complexity. The underlying principle of finite elements is the locality 
of the bases which yields sparse linear systems and offers the possibility of linear 
computational complexity with respect to the number of degrees of freedom. Let 
{Xj I 7 = 1,2,..., A// = dimV//} be such a local basis of Vh- 

We shall derive a basis of the test space Wh defined in (5) by mapping the trial 
basis onto a test basis via some bijective operator a so-called trial-to-test opera¬ 
tor. Due to Assumption (7) such an operator exists, but there are many choices and 
we have to make a design decision. Our choice is that 


Iro = id 


(9) 


which is consistent with almost all existing practical realisations of the method but 
one might as well consider distance minimisation \\{\ — 3X)vh\\v = Ik// ~ 

Wh\\v- 

The condition (9) fixes the (macroscopic) finite element part of 

^vr while the fine scale remainder (1 —Ir)^vr is determined by the variational 
condition in the definition ofWR. Given vr gVr, {I — Ir)^vr e KqyIr satisfies 

Vz G KqyIr : a{z, (I-Ir)^vr) = -a{z,VR). (10) 


This problem is referred to as the fine scale corrector problem for vr eVr. Note 
that Vr can be replaced with any v G V so that (1 —7//)^ can be understood as an 
operator from V into KqyIr. We usually denote this operator the fine scale correction 
operator and write {1 —Ir)^. This operator is the Galerkin projection from 
Vr (or V) into KqyIr related to the adjoint of the sesquilinear form a. It depends on 
the underlying variational problem and equips test functions with problem related 
features that are not present in Vr. In the context of elliptic PDEs, ^ is called the 
finescale Green’s operator lHug95, HS07]. 

For this construction to work we need to assume the well-posedness of the cor¬ 
rector problem (10), i.e., there is some constant j8 > 0 such that 


inf 


AXAi sup I I I I I I I I ^ R ^ AAAA 

07^v€Ker/Ho^w6Ker/H lll^lk |KI|v o^veKer/^ l|v||v ||h'||v' 


a(v, w) 


> J3 < inf 


sup 


a{v,w) 


(11) 
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As for the Galerkin projection Gh onto Vh, these inf-sup conditions do not follow 
from their continuous counterparts (2) (unless a is coercive) and they might hold 
for sufficiently small H only. However, we were able to show in the context of the 
Helmholtz model problem of Section 5 that (11) holds in a much larger regime of the 
discretization parameter H than the corresponding conditions for the standard FEM 
do. In any case, condition (11) implies that the trial-to-test operator ^ = 1 is a 
bounded linear projection operator from V to Wh with operator norm 

ll'^llif(V) = I|1 “ ^\\^(V) = \\^\\^(V) < 

where Ca denotes the continuity constant of the sesquilinear form a. Moreover, 
^\vh • is invertible with = Ir and | j = 1,2,... .Nr} 

defines a basis of Wr with 

1 < j <Nh- 

In general, it cannot be expected (apart from one-dimensional exceptions where 
Ker/// is a broken Sobolev space [HS07]) that the have local support. On 
the contrary, their support will usually be global. However, we will show in the 
next section that they decay very fast in relevant applications; for illustrations see 
Section 4. 

An important special case of the model problem (1) is the hermitian case. Note 
that hermiticity is preserved by the Petrov-Galerkin method in the following sense. 
For any ur^vr ^Vr, it holds that 

a{uH, ^vh) = UH, vh) = a{^ vh, ^uh) = a{vH, ^uh). 

However, this hermiticity is typically lost once ^ is replaced with some approxi¬ 
mation In order to avoid a lack of hermiticity, previous papers such as [MP14] 
mainly used a variant of the method with Wr as the test and trial space. If her¬ 
miticity is important, one should follow this line. In this paper, we trade hermiticity 
for a cheaper method that avoids any costly communication between the fine-scale 
correctors that would be necessary in the hermitian version. 

If the problem is non-hermitian, one might still consider a modified trial space 
based on the adjoint of ^ to improve approximation properties; see [LM09, Mall 1, 
Pet 14] for details. In a setting with a modified trial space, further generalisations 
are possible. Since Vr does not appear any more in the method, its conformity can 
be relaxed as it was recently proposed in [Owhl5] in the context of a multilevel 
solver for Poisson-type problems with LT coefficients. This approach enables one 
to compute very general quantities of the solution such as piecewise mean values. 
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3 Exponential decay of fine-scale correctors 


In many cases, the fine-scale correctors (i.e. the solutions of the fine-scale corrector 
problems (10)) have decay properties better than those of the Green’s function as¬ 
sociated with the underlying full-scale partial differential operator. To elaborate on 
this, we shall now assume that the space V is a closed subspace of with a 

local norm (the notation || • ||v,co means that the V-norm is restricted to some sub- 
domain (0 C ^2). Moreover, the sesquilinear form a is assumed to be local. This is 
the natural setting for scalar second order PDEs. The subsequent arguments can be 
easily generalised to vector-valued problems. 

To be more precise regarding the locality of the basis mentioned above, we shall 
associate the basis functions of Vh with a set of geometric entities called nodes 
(e.g. the vertices of a triangulation) and assume that these nodes are well distributed 
in the domain Q in the sense of local quasi-uniformity. In this context, H refers to 
the maximal distance between nearest neighbours (the mesh size). Given some node 
z G and the corresponding basis function G V//, set the corrector 
and recall from (10) that 

a{w,<j)^) = -a{w,X^) 

for all w G Ker///. 

We aim to show that there are constants c > 0 and C > 0 independent of H and 
R such that 


\\'^^z\\v,n\BR{z) — \\^z\\v,n\BR{z) < ll'^^zllv, (12) 

where (z) denotes the ball of radius R > 0 centred at z. 

We shall show how this result can be established and what kind of assumptions 
have to be made. Let R > 2H and r := R — H > H and let rj G [0,1]) be 

some cut-off function with rj = 0 in Q \ Br{z) , R = lin Br{z), and 

ML-(n)<CnH-^ (13) 

for some generic constant Crf . In general, the fine-scale space KqyIh is not closed 
under multiplication by a cut-off function and we will need to project the truncated 
function back into KqyIh by the operator I— Ir- We assume that the concate¬ 
nation of multiplication by rj and {I — Ir) is stable and quasi-local in the sense 
that 

Vw G KcyIr : 11(1 < Cr]^l^\\w\\y^Bj^i(z)\B^^^^ (14) 

holds with r' \= r — mH and R' := 7? + mH and generic constants > 0 and 
m G No independent of H and z. Although the multiplication by rj is not a stable op¬ 
eration in the full space V (think of a constant function), this result is possible in the 
space of fine scales for example if Ir enjoys quasi-local stability and approximation 
properties; see Section 4 below for an example. The quasi-locality of Ir is also used 
in the next argument. 
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Assuming that the inf-sup condition (11) holds, the corrector satisfies 

< ll(l--fe)((l-»?)<^>z)l|v' 

< /3"'a(w,(l-/H)((l-77)<^>z)) 

= {a{w,(l)^) -a{w,{l -lH){ri(l)^))) 

for some w G Ker/// with ||w||v = 1. Since supp((l —///)((! — C Q \Br{z) 

there is a good chance to actually find a function w with 

suppw c supp((l -/i/)((l c n\Br{z). 


Of course, this is an assumption that needs to be verified in the particular application. 
Under this condition, the term a{w, 0^) = ( 2 (w, A^) vanishes because the supports of 
w and have no overlap. This and (14) imply 




= P- 


'CaQ 


Uh 



\2 


-\m 


2 

v,n\Bj,f{z) 



where Q denotes the continuity constant of the sesquilinear form a. Hence, the 
contraction 

2 ^'2 

UzWv,n\Bi,,{z) ^ 

holds with C' := (jS^^C^Cr^,/^)^. The iterative application of this estimate with R' 

R' — (2m plus relabelling R' R leads to the conjectured decay result (12) 


with constants C := (p^) and c := |log(p^) 


c' 


( 1 ) 

2(2m+l) 


> 0 . 


The exponential decay motivates and justifies the localisation of the fine-scale 
corrector problems to local subdomains of diameter iH where ^ G N is a new dis¬ 
cretization parameter, the so-called oversampling parameter. It controls the pertur¬ 
bation with respect to the ideal global correctors. We will explain this localisation 
procedure on the basis of an example in Section 4 below. As a rule of thumb, 
the localisation to subdomains of diameter iH will introduce an error of order 
^(exp(—^)). As long as this error is small when compared with the inf-sup con¬ 
stant a\\IH\\^^y^^ of the ideal method, the stability and approximation properties of 
the method will be largely preserved. 


4 Application to numerical homogenization of elliptic PDEs 

The first prototypical model problem concerns the diffusion problem — divAVu = f 
in some bounded domain Q cW^ with homogeneous Dirichlet boundary condi¬ 
tion. The difficulty is the strongly heterogeneous and highly varying (non-periodic) 
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diffusion coefficient A. The heterogeneities and oscillations of the coefficient may 
appear on several non-separable scales. We assume that the diffusion matrix A e 
is symmetric and uniformly elliptic with 


0 < a = ess inf inf 


V-V 


Given f e L?'{Q),wq wish to find the unique weak solution u gV := Hq{Q) such 
that 

a{u,v) := [ {AVu)-Vv= [ fv=:F{v) for all v G V. (15) 

Jn Jn 

It is well known that classical polynomial based FEMs can perform arbitrarily badly 
for such problems, see e.g. [BOOO]. This is due to the fact that finite elements tend to 
average unresolved scales of the coefficient and the theory of homogenization shows 
that this way of averaging does not lead to meaningful macroscopic approximations. 
This is illustrated in the introduction. In the simple periodic example of Fig. 1, the 
averaging of the inverse of the diffusion coefficient A (harmonic averaging) would 
have lead to the correct macroscopic representation. 

In computational homogenization, the impact of unresolved microstructures en¬ 
coded in the rough coefficient A on the overall process is taken into account by the 
solution of local microscopic cell problems. While many approaches are empiri¬ 
cally successful and robust for certain multiscale problems, the question whether 
such methods are stable and accurate beyond the strong assumptions of analytical 
homogenization regarding scale separation or even periodicity remained open for 
a long time. Only recently, the existence of an optimal approximation of the low- 
regularity solution space by some arbitrarily coarse generalised finite element space 
(that represents the homogenised problem) was shown in [BLll] and [GGS12]. 
However, the constructions therein include prohibitively expensive global solutions 
of the full fine scale problem or the solution of more involved eigenvalue problems. 
The first efficient and feasible construction, solely based on the solution of localised 
microscopic cell problems, was given and rigorously justified in [MP14] and later 
optimised and generalised in [HP13, HMP15]. A different approach with presum¬ 
ably similar properties was later suggested by [OZB13] along with the notion of 
sparse super-localisation that refiects the locality of the discrete homogenised oper¬ 
ator (similar to the sparsity of standard finite element matrices). 

We shall now explain how the abstract theory of the previous sections is related to 
the method of [MP14] and its variants. Let denote some regular (in the sense of 
Ciarlet) finite element mesh into closed simplices and let Vh := Pi {^h) H V denote 
the space of continuous functions that are affine when restricted to any element 
T e Let Ih ’ V ^ Vh he Si quasi-interpolation operator that acts as a stable 
quasi-local projection in the sense that Ir^Ih = Ih and that for any T ^^r and all 
V G y there holds 


H '||v-///v||^2(7-) + \\Ihv\\v,t < C/^||Vv||v,i2r 


( 16 ) 
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Fig. 3 Standard nodal basis function with respect to the coarse mesh (top left), correspond¬ 
ing ideal corrector 0^ = (top right), and corresponding test basis function = (1 + 
(bottom left). The bottom right figure shows a top view on the modulus of test basis function 
= {l-\-^)Xz with logarithmic color scale to illustrate the exponential decay property. The 
underlying rough diffusion coefficient A is depicted in Fig. 6. 


where Qj refers to some neighbourhood of T (typically the union of T and the ad¬ 
jacent elements) and || • ||y := ||V • One possible choice (among many others) 

is to define 1h := Eh o TIh, where Tin is the piecewise projection onto {^h) 
and Eh is the averaging operator that maps P\ {^h) to Vh by assigning to each inte¬ 
rior vertex the arithmetic mean of the corresponding function values of the adjacent 
elements, that is, for any v G Pi {Wh) and any free vertex z G ^h, 

c^rd{Ke^H :zeK} 


(EHivm = 
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Fig. 4 Element patches Qj/ 
for I = 1,2,3 (from left to 
right) as they are used in the 
localised corrector problem 
(18). 



For this choice, the proof of (16) follows from combining the well-established ap¬ 
proximation and stability properties of Uh and Eh, see for example [DEI2]. The 
choice of 1h in [MP14, HP13] was slightly different. Therein, the -orthogonal 
projection onto Vh played the role of ///. Since this a non-local operator, the analysis 
was based on the fact that the local quasi-interpolation operator of [CV99, Section 
6] has the same kernel and, hence, induces the same method. 

Following the recipe of Section 2.1 and taking into account the present setting 
with an inner product a, the ideal test space Wh '•= (Ker///)^"* is simply the orthog¬ 
onal complement (w.r.t. a) of the fine scale functions Ker///. 

Given the nodal basis of Vh, sl basis of Wh is computed by means of the trial-to- 
test operator ^ = 1 + where 

Vw G Kerin : = —a{A.z^w). (17) 

It is easily checked that the assumptions made in Section 3 are satisfied in the present 
setting. In particular, formula (14) holds with Cr^,/^ = C/^(C/^Cr^ + 1) and m = 2. 
This follows from the product rule, (13), and the local approximation and stability 
properties (16) of ///. This implies the exponential decay as it is stated in (12) with 
constants C and c independent of variations of the diffusion coefficient A. An exam¬ 
ple of a corrector and a test basis function are depicted in Fig. 3 to demonstrate the 
exponential decay. 

We truncate the computational domain of the corrector problems to local subdo¬ 
mains of diameter IH roughly. We have not yet described how to do this in practice. 
The obvious way would be to simply replace in (17) with suitable neighbour¬ 
hoods of the nodes z. This procedure was used in [MP14]. However, it turned out 
that it is advantageous to consider the following slightly more involved technique 
based on element correctors [HP13, HMP15]. 

We assign to any T G its ^-th order element patch Qj.l for a positive integer 
see Fig. 4 for an illustration. Moreover, we define for all v,w G V and CO C Q the 
localised bilinear forms 


j(v,w) := f (AVv)'Vw. 

Jco 


Given any nodal basis function G V//, let (l)z/j ^ KerIn r\WQ'^{QT,i) solve the 
subscale corrector problem 

= -aT{K,w) for all w e Ker/n nWg’^(i2T,e). (18) 
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Fig. 5 Localised element correctors for ^ = 2 and all four elements T adjacent to the vertex 
z = [0.5,0.5] (top), localised nodal corrector 0^^ = = 11t3z^z/,t (bottom left) and corre¬ 
sponding test basis function ^ = (1 (bottom right). The underlying rough dif¬ 

fusion coefficient is depicted in Fig. 6. The computations have been performed by standard linear 
finite elements on local fine meshes of with h = 2~^. See Fig. 3 for a comparison with the ideal 
global corrector and basis. 


Let (l)z^£ := and define the test function 

^Z-I^ • — ^z • 

The localised test basis function and the underlying correctors can be 
seen in Fig. 5. Note that we impose homogeneous Dirichlet boundary condition on 
the artificial boundary of the patch which is well justified by the fast decay. 

More generally, we may define the localised correction operator by 

'•= ^ ^h{z)(I)z,£ 
ze^H 

as well as the localised trial-to-test operator 

^£Vh := 1 -\-'^£Vh = ^ vh{z)Az^£. 

ze^H 


The space of test functions then reads 
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Wh := ^(Vh = span{Aj^<> : z € 

and the (localised) multiscale Petrov-Galerkin FEM seeks C Vh such that 

= {f for all wh/ C (19) 

In previous papers [MP14, HP13, HMP15] we have considered the symmetric ver¬ 
sion with Wh^i as trial and test space and also the reverse version with as the 
trial space and Vh as test space [EGH15]. All these methods are essentially equal in 
the ideal case and there are no major changes in the output after localisation (when 
only the Vh part of the discrete solution is considered). When it comes to imple¬ 
mentation and computational complexity, the present Petrov-Galerkin version has 
the advantage that there is no communication between the correctors. This means 
that the fine-scale solutions of the corrector problems need not to be stored but only 
their interaction with the standard nodal basis functions in their patches; see 

also [EGH15] for further discussions regarding those technical details. 

The error analysis of the localised method follows similar arguments. Let uh ^ 
Vh be the ideal Petrov-Galerkin approximation and let en '=uh — uh/ G Vh denote 
the error with respect to the ideal method. Then there exists some zh ^ with 
||z//||v = 1 such that 

(X 

TT^-T- W^hWv < Cl{eH,ZH) = Cl{uH^i-U,ZH-ZH/)i 

¥h\\^(v) 

where zh,i ^ The exponential decay property allows one to choose zh/ in 
such a way that \\zh—Zh/\\v < Cexp(—c^); see for instance [HP13, HMP15]. This 
shows that 

\\u — Uh^i\\v < \\u — Uh\\v-V\\uh—Uh/\\v 

< \\u — Uh\\v + - CQ^XZ^{—ct)\\u — UH^l\\V‘ 

We shall emphasise that, in the present context, the constants C and c are inde¬ 
pendent of variations of the rough diffusion tensor but they may depend on the 
contrast (the ratio between the global upper and lower bound of A). Using (8), this 
shows that the moderate choice I > \ \og{a/{2\\lH\\^{y)^aC))\/c = ^(1) implies 
the quasi-optimality (and also the well-posedness) of the Petrov-Galerkin method 
with respect to the V-norm 

\\u-uh/\\v <AVh\\^(v) min \\u-vh\\v^ 

vh^Vh 

With regard to the fact that the V-best approximation may be poor and standard 
Galerkin would have provided us with an even better estimate at lower cost, this 
result is maybe not very impressive. Let us see if we can do something similar for the 
L^-norm which appears to be the relevant measure in the context of homogenization 
problems. A standard duality argument shows that 
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Fig. 6 Diffusion coefficient in the numerical experiment of Section 4 and coarsest mesh. 


= a{eH,ZH) =a{u-UH,e,ZH-ZH,e) 

for some zh € Wh with ||zh||v < Zh/ ■= ^(^hZh € 

Wh/- Similar arguments as before yield 

l|M-MH.t|lz, 2 (^) <Ci min \\u-Vh\\l 2 (^q)+C 2 Cw{-c£) min \\u-Vh\\v, 

vH ^ ^ H ^ 

where C\ := 114^11^(^2(12)) ^2 •= CaCC^OC~^\\lH\\^(v)' This shows that the 

method is accurate also in the L^-norm regardless of the regularity of u. If the over- 
sampling parameter is chosen such that I > \ogH, then the method is (^{H) accurate 
in L^{Q) with no pre-asymptotic phenomena. This is the best worst-case rate one 
can expect for general / G V' and A e L°°. 

Note that the previous results hold true for general L°°-coefficients and all con¬ 
stants are independent of the variations of the diffusion tensor as far as the contrast 
remains moderately bounded. In particular, the approach is by no means restricted 
to periodic coefficients or scale separation. For a more detailed discussion of high- 
contrast problems in this context we refer to [PS 15]. 

The final step towards a fully practical method is the discretization of the fine- 
scale corrector problems. With regard to the possible low regularity of the solution, 
Pi finite elements on a refined mesh appears reasonable, but any other type of 
discretisation is possible. Obviously, the fine-scale discretization parameter h has to 
be chosen fine enough to resolve all relevant features of the diffusion coefficient. 
The previous theory can be transferred to this case in a straight-forward way and we 
refer to [MP14, HP13, HM14] for the technical details. 

To illustrate the previous estimates, we close this section with a numerical exper¬ 
iment. Let Q be the unit square and the outer force f=linQ. Consider the coeffi¬ 
cient A that is piecewise constant with respect to a uniform Cartesian grid of width 
2~^. Its values are randomly chosen between 1 and 10; see Fig. 6. Consider uniform 
coarse meshes of size Ff = 2“^,..., 2~^ of Q that certainly do not resolve 

the rough coefficient A appropriately. The reference mesh % has width h = 2~^. 
Since no analytical solutions are available, the standard finite element approxima¬ 
tion Uh G Vh on the reference mesh % serves as the reference solution. Doing this. 









Variational Multiscale Stabilization and the Exponential Decay of Fine-scale Correctors 


17 


Fig. 7 Numerical experiment 
of Section 4. Relative L^- 
errors of multiscale Petrov- 
Galerkin FEM (19) versus 
the number of degrees of 
freedom Nh ~ where 
// = , 2~^ is the uni¬ 

form coarse mesh size. The 
localisation parameter varies 
between ^ = 1,..., 3. The 
Pi -FE solution and the best- 
approximation in the Pi-FE 
space on the same coarse 
meshes are depicted for com¬ 
parison. 


10 ° 10 ^ 10 ^ 10 ° 
Nh»H-2 



we assume that Uh is sufficiently accurate and, necessarily, that % resolves the dis¬ 
continuities of A. The corrector problems are also are also solved on this scale of 
numerical resolution. 

The numerical results, i.e. errors with respect to the reference solution Uh are 
depicted in Fig. 7. The results are in agreement with the theoretical results. They 
are even better in the sense that £ = 1 seems to be sufficient for quasi-optimality 
(with respect to Uh) in the present setup and parameter regime. We expect that the 
true errors with respect to u would behave similar in the beginning but level off at 
some point when the reference error starts to dominate the upscaling error. Still, 
the experiment clearly indicates that numerical homogenization is possible for very 
general -coefficients. 

We refer to [MP14, HP13, HMP15, HMP14b, EGH15, AH14, BP14, EGMP13, 
HM14, MP15b, HMP14a] for many more numerical experiments for several model 
problems including nonlinear stationary and non-stationary problems as well as 
eigenvalue problems. 


5 Application to high-frequency acoustic scattering 

This section will show that the abstract framework of Section 2-3 is indeed applica¬ 
ble beyond the coercive and symmetric model problem of the previous section. We 
consider the scattering of acoustic waves at a sound-soft scatterer modelled by the 
Helmholtz equation over a bounded Lipschitz domain C (J = 1,2,3), 

—Au—K^u = f in^2. 


(20.a) 
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along with mixed boundary conditions of Dirichlet and Robin type 


u = 0 on I]), (20.b) 

'Vu’V — iKu = 0 on Fr. (20.c) 

Here, the wave number /c ^ 1 is real and positive, i denotes the imaginary unit and 
f e L^{Q^C). We assume that the boundary F := dQ consists of two components 


dQ=FDUFR^ F/)nF/^ = 0 


where Fq encloses the scatterer and Fr is an artificial truncation of the whole un¬ 
bounded space. The vector v denotes the unit normal vector that is outgoing from 

Q. 

Given / G L^(^2,C), we wish to find w G V := {v G | v = 0 on F^} 

such that, for all v G V, 

a{u^v) := / Vw-Vv— K*^/ uv — iK uv = / fv=:F{w). (21) 

Jn Jn Jtr Jn 

The space V is equipped with the usual /c-weighted norm ||v||y := + 

||Vv||^ 2 (^)- The presence of the Robin boundary condition (20.c) ensures that this 
variational problem is well-posed in the sense of (2) with a = l/Cst(K‘) for some 
K*-dependent stability constant Cst(K:); see e.g. [EM12]. The dependence on the 
wave number K is not known in general. An exponential growth with respect to 
the wave number is possible [BCWG+11] in non-generic domains. In most cases, 
the growth seems to be only polynomially, although this is an empirical rather than 
a theoretical statement, and sufficient geometric conditions for this to hold are rare 
[EM 12, Mel95, CE06, MIB96]. Eor the above scattering problem, we know that 
< ^{k) if Q is convex and if the scatterer is star-shaped [Het07]. 

It is this K:-dependence in the stability of the problem that makes the numeri¬ 
cal approximation by EEM or related schemes extremely difficult in the regime of 
large wave numbers. Any perturbation of the problem, e.g. by some discretization, 
can be amplified by Cst(K*). We have seen in the introduction that this is indeed 
observed in practice and causes a pre-asymptotic effect known as the pollution ef¬ 
fect or numerical dispersion [BSOO]. This effect puts very restrictive assumptions 
on the smallness of the underlying mesh that is much stronger than the minimal 
requirement for a meaningful representation of highly oscillatory functions from 
approximation theory, that is, to have at least 5 — 10 degrees of freedom per wave 
length and coordinate direction. 

It is the aim of many newly developed methods to overcome or at least to re¬ 
duce the pollution effect; see [TE06, EW09, EWll, HMPll, ZMD+11, DGMZ12] 
among many others. However, the only theoretical results regard high-order EEMs 
with the polynomial degree p coupled to the wave number K via the relation 
p ^ logK [MSIO, MSll, MPS13, EM12]. Under this moderate assumption, those 
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methods are stable and quasi-optimal in the regime HKjp < 1 for certain model 
Helmholtz problems. 

The multiscale method of [Pet 14] then showed that pollution in the numerical 
approximation of the Helmholtz problem can also be cured for a fairly large class of 
Helmholtz problems, including the acoustic scattering from convex non-smooth ob¬ 
jects, by stabilization in the present framework. If the data of the problem (domain, 
boundary condition, force term) allows for polynomial-in-/c bounds of Cst(K‘) and if 
the resolution condition Hk <1 and the oversampling condition log(K:)/^ < 1 are 
satisfied, then the method is stable and quasi-optimal in the V-norm. 

The recent paper [GP15] interprets the method of [Pet 14] in the present frame¬ 
work and we recall it here very briefiy. Given the same discrete setup as in the 
previous section with some simplicial mesh corresponding Pi FE space Vh '= 
P\{^h) nV, and quasi-interpolation operator Ir - V the multiscale Petrov- 

Galerkin method is formally defined in the same way. We simply replace the inner 
product of Section 4 with the sesquilinear form a of this section. 

Given any nodal basis function G V//, we construct a corresponding test basis 
function by the same procedure as in the previous section, ;= + 0^^^, 

where := and solves the cell problem 

= —ciriw, A^) for all w G KerI r with suppw C Qt- 

Here, 

•= / Vt/-Vv—K*^/ uv — iK uv 

jQr\(o JQr\(o JrRr\d(o 

for (O G T }. Note that the corrector problem inherits the boundary condition 

from the original problem when the patch boundary coincides with the boundary 
of . On the part of the patch boundary that falls in the interior of , we simply 
put the homogeneous Dirichlet condition. A major observation is that this corrector 
problem is well-posed and, in particular, coercive with j8 = 1/3 under the condition 
Hk< Cres for some given constant 0 < Cres = ^(1) that only depends on the constant 
in (16) but not on H or K. This is because a satisfies a Garding inequality and fine- 
scale functions satisfy ||w||^ 2 (^) < This coercivity also implies 

the desired exponential decay of the ideal correctors so that the choice Qt/ is well 
justified. This can also be observed in Fig. 8. 

The space of localised test functions then reads Wr^£ := span{A^^^ : z G ^r} and 
the multiscale Petrov-Galerkin FEM seeks ur^£ G Vr such that 

(^{ur^£^wr^£) = F{wr^£) for all wr^£ G Wr^£. (22) 

The quasi-optimality result of the previous section is easily transferred to the 
present setting. The resolution condition Hk < Cres and the oversampling condition 

£> |log(a/(2||/H||^(v)QC))|/c= ^(logQ,(jc)) 

imply the quasi-optimality (and stability) of the multiscale Petrov-Galerkin method 
with respect to the V-norm 
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Fig. 8 Real and imaginary part of the ideal corrector (top left and middle). The top right 
figure shows a top view on the modulus of test basis function = with logarithmic 

color scale to illustrate the exponential decay property. The underlying computational domain is 
the unit square with a Robin boundary condition everywhere. The wave number K" = 2^ is chosen 
such that the resolution condition on the coarse mesh is just satisfied. The localised nodal corrector 
(bottom left) and corresponding test basis function ^ = £XXz (bottom right) are real¬ 
valued because the patch boundary doesn’t touch the domain boundary. The local fine meshes used 
in the computation have width h = 2~^. 


< 2||///||^(y) min \\u — vh\\v‘ 
vh^Vh 

Here, the constants c and C are related to the exponential decay of the test basis 
(cf. (12)) and they are independent of K under the resolution condition. We shall 
emphasise that such a best-approximation property does not hold for standard FEMs 
which require e.g. k^H < 1 for quasi-optimality [Mel95] in the case of pure Robin 
boundary conditions on a convex planar domain. The FEM approximation is not 
even known to exist unless < 1 in the simplest model problem without a 

scatterer [Wul4]. 

For the multiscale Petrov-Galerkin method, the result means that pollution effects 
do not occur. Note that the resolution condition Hk < Cres is somewhat minimal. 
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Fig. 9 Computational domain 
of the model problem of 
Section 5 and coarsest mesh. 



because any meaningful approximation of the highly oscillatory solution of (20) 
requires at least 5 — 10 degrees of freedom per wave length and coordinate direction. 
Saying this, we assume that the fine scale corrector problems are solved sufficiently 
accurate; see [GP15, Pet 14] for details. 

We shall present a numerical experiment taken from [Pet 14] where this ver¬ 
sion of the method was already considered experimentally. Consider the scatter¬ 
ing from sound-soft scatterer occupying the triangle Qd- The Sommerfeld radia¬ 
tion condition of the scattered wave is approximated by the Robin boundary con¬ 
dition on the boundary Fr := dQ.R of the unit square so that := (0,1)^ \ is 
the computational domain; see Fig. 9. Given the wave number k = 2^, the incident 
wave Uinc{x) := exp(//c v- [cos(0.5), sin(0.5)]^) is prescribed via an inhomogeneous 
Dirichlet boundary condition on Fd := BQd and the scattered wave satisfies (20.a) 
with / = 0 and the boundary conditions 

u = on Tj) , 

Vu' V — iKu = 0 on Fr. 

The error analysis extends to this setting in a straight-forward way. 

We choose uniform coarse meshes of widths H = 2“^,... ,2“^ as depicted in 
Fig. 9. The reference mesh % is derived by uniform mesh refinement of the coarse 
meshes and has mesh width h = 2~^. The corresponding Pi conforming finite ele¬ 
ment approximation on the reference mesh is denoted by Vh. As in the previous 
section, we compare the coarse scale approximations G Vh with some refer¬ 
ence solution Uh GVh. 

Fig. 10 depicts the results for the multiscale Petrov-Galerkin method and shows 
that the pollution effect that is present in the Pi FEM is eliminated when £ is mod¬ 
erately increased. For the present wave number ^ = 2 is sufficient. 

Further numerical experiments are reported in [Pet 14] and [GP15]. It is worth 
noting that the latter work also exploits the homogeneous structure of the PDF co¬ 
efficients in the sense that only very few of the fine-scale corrector problems are ac¬ 
tually solved due to translation invariance and symmetry. This makes the approach 
competitive. 

A very natural and straight forward generalisation of the method would be the 
case of heterogeneous media. The previous section plus the analysis of this section 
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Fig. 10 Numerical experi¬ 
ment of Section 5: Relative 
V -norm errors of multiscale 
Petrov-Galerkin method (22) 
with wave number k = 2^ de¬ 
pending on the number of de¬ 
grees of freedom Nh ~ 
where H = 2 ~^,..., 2~^ is 
the uniform coarse mesh 
size. The reference mesh size 
h = 2~^ remains fixed. The 
oversampling parameter £ 
varies between 1 and 3. The 
Pi -FE solution and the best- 
approximation in the Pi-FE 
space on the same coarse 
meshes are depicted for com¬ 
parison. 



Strongly indicate the potential of the method to treat high oscillations or jumps in 
the PDE coefficients and the pollution effect in one stroke. 


6 Final remarks 

We have presented an abstract framework for the stabilization of numerical meth¬ 
ods for multiscale partial differential equations with some focus on highly oscil¬ 
latory problems. The methodology is based on the variational multiscale method 
and the more recent development of localised orthogonal decompositions. We have 
provided an abstract numerical analysis of the method which is applied to two rep¬ 
resentative model problems, a homogenization problem and a scattering problem. 
We have shown that the methodology can indeed eliminate critical scale-dependent 
pre-asymptotic effects in these cases. While the framework has already been applied 
successfully to other problem classes such as linear and non-linear eigenvalue prob¬ 
lems, we expect that the framework will also be useful for convection-dominated 
fiow, the problem that the variational method was initially designed for. 

The multiscale method presented in this paper is shown to be stable and accurate 
under moderate assumptions on the discretization parameters relative to character¬ 
istic parameters and length scales of the problem. These valuable properties require 
the pre-computation of the test basis on subgrids. These pre-computations are both 
local and independent, but the worst-case (serial) complexity of the method can ex¬ 
ceed the cost of a direct numerical simulation on a global sufficiently fine mesh. 
If the inherent parallelism of the local cell problems cannot be exploited during 
the computation, we still expect a significant gain with respect to computational 
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complexity if the pre-computation can be reused several times in the context of pa¬ 
rameter studies, coupled problems, optimal control problems or inverse problems. 
In many cases, there is also a lot of redundancy in the local problems which allows 
one to reduce the number of local problems drastically as it is shown in [GP15] in 
the context of acoustic scattering. We expect that this technique can be generalised 
to far more general situations using modern techniques of model order reduction 
[RHP08, AB14]. 

We may close the discussion with some rather philosophical remark regarding 
the stabilization of FEMs and their inter-element continuity properties. Presently, it 
is believed, e.g. in the context of time-harmonic wave propagation, that stability can 
be increased by relaxing inter-element continuity within a discontinuous Galerkin 
(DG) framework. The large number of variants includes the ultra weak variational 
formulation [CD98], Trefftz methods [HMPll], DPG [ZMD+11, DGMZ12], or 
the continuous interior penalty method [Wul4]. There may be some truth in this 
but the general impression that relaxing continuity is the only way is certainly 
false as one can observe from the method presented in this paper. The multiscale 
Petrov-Galerkin does quite the opposite. The regularity of the test functions is in¬ 
creased compared to standard continuous finite elements, because they are solutions 
of second-order elliptic problems (at least in the ideal case). In general, test func¬ 
tions wh G Wh have the property that divAVw// G In the context of the 

Helmholtz model problem of Section 5 where A = 1 this means that AWh C 
If is convex and boundary conditions are appropriate, then Wh C (this 

can be observed for one basis function in Fig. 8). In this respect, our methodol¬ 
ogy clearly indicates that increased differentiability might as well lead to increased 
stability and accuracy. Similar effects have been observed for eigenvalue compu¬ 
tations in IGA [CRBH06] and also LOD [MP15b]. This shows that breaking the 
inter-element continuity is not at all necessary for stability. 
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